home *** CD-ROM | disk | FTP | other *** search
/ CD-ROM Today - The Disc! 5 / CD-ROM Today - The Disc (Issue 5)(November 1994).ISO / mac / Mac shareware / Education / RLaB / test < prev    next >
Text File  |  1994-09-21  |  47KB  |  1,822 lines
  1. tic();
  2. //
  3. //  RLaB test program
  4. //
  5.  
  6. //
  7. // Test Parameters and some functions we will need later
  8. //
  9.  
  10. X = 3;
  11.  
  12. epsilon = function() 
  13. {
  14.   eps = 1.0;
  15.   while((eps + 1.0) != 1.0) 
  16.   {
  17.     eps = eps/2.0;
  18.   }
  19.   return eps;
  20. };
  21.  
  22. eps = epsilon();
  23.  
  24. eye = function( m , n ) 
  25. {
  26.   local(i, N, new);
  27.  
  28.   if (!exist (n))
  29.   {
  30.     if(m.n != 2) { error("only 2-el MATRIX allowed as eye() arg"); }
  31.     new = zeros (m[1], m[2]);
  32.     N = min ([m[1], m[2]]);
  33.   else
  34.     if (class (m) == "string" || class (n) == "string") {
  35.       error ("eye(), string arguments not allowed");
  36.     }
  37.     if (max (size (m)) == 1 && max (size (n)) == 1)
  38.     {
  39.       new = zeros (m[1], n[1]);
  40.       N = min ([m[1], n[1]]);
  41.     else
  42.       error ("matrix arguments to eye() must be 1x1");
  43.     }
  44.   }
  45.   for(i in 1:N)
  46.   {
  47.     new[i;i] = 1.0;
  48.   }
  49.   return new;
  50. };
  51.  
  52. symm = function( A )
  53. {
  54.   return (A + A')./2;
  55. };
  56.  
  57. //-------------------------------------------------------------------//
  58.  
  59. // Synopsis:    Pascal matrix.
  60.  
  61. // Syntax:    P = pascal ( N )
  62.  
  63. // Description:
  64.  
  65. //    The Pascal matrix of order N: a symmetric positive definite
  66. //    matrix with integer entries taken from Pascal's triangle. The
  67. //    Pascal matrix is totally positive and its inverse has integer
  68. //    entries.  Its eigenvalues occur in reciprocal pairs. COND(P)
  69. //    is approximately 16^N/(N*PI) for large N. PASCAL(N,1) is the
  70. //    lower triangular Cholesky factor (up to signs of columns) of
  71. //    the Pascal matrix.   It is involutary (is its own
  72. //    inverse). PASCAL(N,2) is a transposed and permuted version of
  73. //    PASCAL(N,1) which is a cube root of the identity.
  74.  
  75. //      References:
  76. //      R. Brawer and M. Pirovino, The linear algebra of the Pascal matrix,
  77. //           Linear Algebra and Appl., 174 (1992), pp. 13-23 (this paper
  78. //           gives a factorization of L = PASCAL(N,1) and a formula for the
  79. //           elements of L^k).
  80. //      S. Karlin, Total Positivity, Volume 1, Stanford University Press,
  81. //           1968.  (Page 137: shows i+j-1 choose j is TP (i,j=0,1,...).
  82. //                   PASCAL(N) is a submatrix of this matrix.)
  83. //      M. Newman and J. Todd, The evaluation of matrix inversion programs,
  84. //           J. Soc. Indust. Appl. Math., 6(4):466--476, 1958.
  85. //      H. Rutishauser, On test matrices, Programmation en Mathematiques
  86. //           Numeriques, Editions Centre Nat. Recherche Sci., Paris, 165,
  87. //           1966, pp. 349-365.  (Gives an integral formula for the
  88. //           elements of PASCAL(N).)
  89. //      J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
  90. //           Birkhauser, Basel, and Academic Press, New York, 1977, p. 172.
  91. //      H.W. Turnbull, The Theory of Determinants, Matrices, and Invariants,
  92. //           Blackie, London and Glasgow, 1929.  (PASCAL(N,2) on page 332.)
  93.  
  94. //    This file is a translation of pascal.m from version 2.0 of
  95. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  96. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  97.  
  98. //-------------------------------------------------------------------//
  99.  
  100. pascal = function ( n , k )
  101. {
  102.   local(P, i, j, k, n)
  103.  
  104.   if (!exist (k)) { k = 0; }
  105.  
  106.   P = diag( (-1).^[0:n-1] );
  107.   P[; 1] = ones(n,1);
  108.  
  109.   //  Generate the Pascal Cholesky factor (up to signs).
  110.  
  111.   for (j in 2:n-1)
  112.   {
  113.     for (i in j+1:n)
  114.     {
  115.       P[i;j] = P[i-1;j] - P[i-1;j-1];
  116.     }
  117.   }
  118.  
  119.   if (k == 0)
  120.   {
  121.     P = P*P';
  122.   else if (k == 2) {
  123.     P = rot90(P,3);
  124.     if (n/2 == round(n/2)) { P = -P; }
  125.   }}
  126.  
  127.   return P;
  128. };
  129.  
  130. //
  131. //-------------------- Test Relational Expressions -------------------
  132. //
  133. printf("\tstart scalar tests...\n");
  134. printf("\tstart relational tests...\n");
  135.  
  136. //    SCALAR CONSTANTS (REAL)
  137. if( !(1<2) ) { error(); }
  138. if( !(1<=2) ) { error(); }
  139. if( 1>2 ) { error(); }
  140. if( 1>=2 ) { error(); }
  141. if( 1==2 ) { error(); }
  142. if( !(1!=2) ) { error(); }
  143. if( !1 ) { error(); }
  144. if( !!!1) { error(); }
  145.  
  146. if( !([1]<[2]) ) { error(); }
  147. if( !([1]<=[2]) ) { error(); }
  148. if( [1]>[2] ) { error(); }
  149. if( [1]>=[2] ) { error(); }
  150. if( [1]==[2] ) { error(); }
  151. if( !([1]!=[2]) ) { error(); }
  152. if( ![1] ) { error(); }
  153. if( !!![1]) { error(); }
  154.  
  155. //    SCALAR CONSTANTS (COMPLEX)
  156. if( !(1+2i<2+3i) ) { error(); }
  157. if( !(1+2i<=2+3i) ) { error(); }
  158. if( 1+2i>2+3i ) { error(); }
  159. if( 1+2i>=2+3i ) { error(); }
  160. if( 1+2i==2+3i ) { error(); }
  161. if( !(1+2i!=2+3i) ) { error(); }
  162. if( !1+2i ) { error(); }
  163. if( !!!1+2i) { error(); }
  164.  
  165. if( !([1+2i]<[2+3i]) ) { error(); }
  166. if( !([1+2i]<=[2+3i]) ) { error(); }
  167. if( [1+2i]>[2+3i] ) { error(); }
  168. if( [1+2i]>=[2+3i] ) { error(); }
  169. if( [1+2i]==[2+3i] ) { error(); }
  170. if( !([1+2i]!=[2+3i]) ) { error(); }
  171. if( ![1+2i] ) { error(); }
  172. if( !!![1+2i]) { error(); }
  173.  
  174. //    SCALAR ENTITIES (REAL)
  175. a=1;b=2;
  176. if( !(a<b) ) { error(); }
  177. if( !(a<=b) ) { error(); }
  178. if( a>b ) { error(); }
  179. if( a>=b ) { error(); }
  180. if( a==b ) { error(); }
  181. if( !(a!=b) ) { error(); }
  182. if( !a ) { error(); }
  183. if( !!!a) { error(); }
  184.  
  185. if( !([a]<[b]) ) { error(); }
  186. if( !([a]<=[b]) ) { error(); }
  187. if( [a]>[b] ) { error(); }
  188. if( [a]>=[b] ) { error(); }
  189. if( [a]==[b] ) { error(); }
  190. if( !([a]!=[b]) ) { error(); }
  191. if( ![a] ) { error(); }
  192. if( !!![a]) { error(); }
  193.  
  194. //    SCALAR ENTITIES (COMPLEX)
  195. a=1+2i;b=2+3i;
  196. if( !(a<b) ) { error(); }
  197. if( !(a<=b) ) { error(); }
  198. if( a>b ) { error(); }
  199. if( a>=b ) { error(); }
  200. if( a==b ) { error(); }
  201. if( !(a!=b) ) { error(); }
  202. if( !a ) { error(); }
  203. if( !!!a) { error(); }
  204.  
  205. if( !([a]<[b]) ) { error(); }
  206. if( !([a]<=[b]) ) { error(); }
  207. if( [a]>[b] ) { error(); }
  208. if( [a]>=[b] ) { error(); }
  209. if( [a]==[b] ) { error(); }
  210. if( !([a]!=[b]) ) { error(); }
  211. if( ![a] ) { error(); }
  212. if( !!![a]) { error(); }
  213.  
  214. if (! all (all (-2 <  rand(4,4)))) { error (); }
  215. if (! all (all (-2 <= rand(4,4)))) { error (); }
  216. if (! all (all ( 2 >  rand(4,4)))) { error (); }
  217. if (! all (all ( 2 >= rand(4,4)))) { error (); }
  218.  
  219. if (! all (all (rand(4,4) >  -2))) { error (); }
  220. if (! all (all (rand(4,4) >= -2))) { error (); }
  221. if (! all (all (rand(4,4) <   2))) { error (); }
  222. if (! all (all (rand(4,4) <=  2))) { error (); }
  223.  
  224. if (! all (all (rand (4,4) >  -rand (4,4)))) { error (); }
  225. if (! all (all (rand (4,4) >= -rand (4,4)))) { error (); }
  226. if (! all (all (-rand (4,4) <  rand (4,4)))) { error (); }
  227. if (! all (all (-rand (4,4) <= rand (4,4)))) { error (); }
  228.  
  229. //------------------------- Test REAL SCALARS ------------------------
  230. //
  231. //    CONSTANTS
  232. //      Addition
  233. if(1+2 != 3) { error(); }
  234. //      Subtraction
  235. if(1-2 != -1) { error(); }
  236. //      Multiply
  237. if(1*2 != 2) { error(); }
  238. //      Divide
  239. if(1/2 != 0.5) { error(); }
  240. //      Power
  241. if(2^2 != 4) { error(); }
  242. if(4^0 != 1) { error(); }
  243. //      Unary Minus
  244. if(2-3 != -1) { error(); }
  245. //
  246. //    ENTITIES
  247. //
  248. a = 1; b = 2; c = 3; d = 0.5;
  249. //      Addition
  250. if(a+b != c) { error(); }
  251. //      Subtraction
  252. if(a-b != -a) { error(); }
  253. //      Multiply
  254. if(a*b != b) { error(); }
  255. //      Divide
  256. if(a/b != d) { error(); }
  257. //      Power
  258. if(b^b != b*b) { error(); }
  259. if(b^0 != 1) { error (); }
  260. //      Unary Minus
  261. if(b-c != -a) { error(); }
  262. //
  263. //  ENTITIES & CONSTANTS
  264. //
  265. if(a+2 != c) { error(); }
  266. //      Subtraction
  267. if(1-b != -a) { error(); }
  268. //      Multiply
  269. if(a*2 != 2) { error(); }
  270. //      Divide
  271. if(1/b != d) { error(); }
  272. //      Power
  273. if(2^b != b*b) { error(); }
  274. //      Unary Minus
  275. if(b-3 != -a) { error(); }
  276. //
  277. //------------------------Test COMPLEX SCALARS -------------------------
  278. //
  279. //    CONSTANTS
  280. if(sqrt(-1) != 1i) { error(); }
  281. //      Addition
  282. if((1+2i)+(2+3i) != (3+5i)) { error(); }
  283. //      Subtraction
  284. if((1+2i)-(3+4i) != (-2-2i)) { error(); }
  285. //      Multiply
  286. if((1+2i)*(3+4i) != (-5+10i)) { error(); }
  287. //      Divide
  288. if((1+2i)/(3-4i) != (-.2+.4i)) { error(); }
  289. //      Power
  290. //      Precision problems prevent us from testing these. Have to
  291. //      be checked by hand.
  292. //  (1+2i)^2 = -3 + 4i
  293. //  (1+2i)^.5 = 1.272 + 7.862e-1i
  294. //  if((1+2i)^2 != (-3+4i)) { error(); }
  295. //  if((1+2i)^10 != (237+3116i)) { error(); }
  296. //      Unary Minus
  297. if(-(1+2i) != -1-2i) { error(); }
  298. //
  299. //    ENTITIES
  300. //
  301. a = 1+2i; b = 3+4i; c = 4+6i;
  302. if(a+b != c) { error(); }
  303. //      Subtraction
  304. if(a-b != -2-2i) { error(); }
  305. //      Multiply
  306. if(a*b != -5+10i) { error(); }
  307. //      Divide
  308. //if(a/(3-4i) != -.2+.4i) { error(); }
  309. //      Power
  310. //  if(b^b != b*b) { error(); }
  311. //      Unary Minus
  312. if(-a != -1-2i) { error(); }
  313. //
  314. //    ENTITIES & CONSTANTS
  315. //
  316. if(a+(3+4i) != c) { error(); }
  317. //      Subtraction
  318. if((1+2i)-b != -2-2i) { error(); }
  319. //      Multiply
  320. if(a*(3+4i) != -5+10i) { error(); }
  321. //      Divide
  322. //if(a/(3-4i) != -.2+.4i) { error(); }
  323. //      Power
  324. //if(b^b != b*b) { error(); }
  325. //      Unary Minus
  326. if(-(1+2i) != -1-2i) { error(); }
  327. //
  328. // String - Numerical Equalities
  329. //
  330. if ((1 == "1")) { error(); }
  331. if (([1] == "1")) { error(); }
  332. if (("1" == 1)) { error(); }
  333. if (("1" == [1])) { error(); }
  334.  
  335. if (rand(3,3) == "str") { error(); }
  336. if ("str" == rand(3,3)) { error(); }
  337. if (!any (any ((["1", "2"; "3", "4"] == "1") == [1,0;0,0]))) { error(); }
  338.  
  339. //
  340. //----------------------- Test REAL MATRICES ---------------------------
  341. //
  342. printf("\tstart matrix tests...\n");
  343. printf("\t\treal-matrices\n");
  344. //  Read in test matrices
  345. //
  346. read("test.input");
  347. //
  348. //  Matrix construction
  349. //
  350. if(all([3] != matrix(3))) { 
  351.   error(); 
  352. }
  353. if(any([1;2;3] != [1,2,3]')) {
  354.   error();
  355. }
  356. if(any (any (m0 != zeros(2,2)))) { error(); }
  357. if(any (any (m1 != 1+zeros(2,2)))) { error(); }
  358. if(any (any (m2 != [1,2;3,4]))) { error(); }
  359. if(any (any (m3 != [1+2i,2+3i;3+4i,5+6i]))) { error(); }
  360. if(any (any ([1,2;3+0i,4+0i] != m2))) { error(); }
  361. if(any (any ([m2] != m2))) { error(); }
  362. //
  363. //  Matrix indexing
  364. //
  365. p = pascal(6);
  366. if (!all (all (p[ [1:3] ; ] == p[ [1:3]' ; ]))) { error (); }
  367. if (!all (all (p[ ; [1:3] ] == p[ ; [1:3]' ]))) { error (); }
  368. if (!all (all (p[ [2:6] ; [2:6] ] == p[ [2:6]'; [2:6]' ]))) { error (); }
  369. if (!all (all (p[ [2:6]' ; [2:6] ] == p[ [2:6]'; [2:6]' ]))) { error (); }
  370. if (!all (all (p[ [6:12]' ] == p[ [6:12] ]))) { error (); }
  371.  
  372. //
  373. //  Sub-Matrix promotion
  374. //
  375. if(m2[2;2] != 4) { error(); }
  376. if(any(m2[2;] != [3,4])) { error(); }
  377. if(any(m2[;2] != [2,4]')) { error(); }
  378. i=2;j=1;
  379. if(m2[i;j] != 3) { error(); }
  380. i=1;j=2;
  381. if(m2[i;j] != 2) { error(); }
  382. m = [1,2,3;4,5,6;7,8,9];
  383. if(any(m[1;1,2] != [1,2])) {
  384.   error();
  385. }
  386. if(any(m[1,2;1] != [1;4])) {
  387.   error();
  388. }
  389. if(any (any (m[1,2;1,2] != [1,2;4,5]))) {
  390.   error();
  391. }
  392. //
  393. if(m3[2;2] != (5+6i)) { error(); }
  394. if(any(m3[2;] != [3+4i,5+6i])) { error(); }
  395. if(any(m3[;2] != conj([2+3i,5+6i]'))) { error(); }
  396. //
  397. //  Automatic creation, extension
  398. //
  399. if(any (any ((mm[3;3]=10) != [0,0,0;0,0,0;0,0,10]))) { error(); }
  400. a=[1,2,3;4,5,6];
  401. if(any (any ((a[3;1]=10) != [1,2,3;4,5,6;10,0,0]))) { error(); }
  402. a=[1,2;3,4];
  403. if(any (any ((a[3,4;3,4]=[5,6;7,8]) != [1,2,0,0;3,4,0,0;0,0,5,6;0,0,7,8]))) {
  404.   error();
  405. }
  406. mmm[2;] = a[1;];
  407. //
  408. //  Matrix binary operations
  409. //
  410. a = m2; b = [5,6;7,8];
  411. if(any (any (a+a != [2,4;6,8]))) { error(); }
  412. if(any (any (a-a != zeros(2,2)))) { error(); }
  413. if(any (any (2+a != [3,4;5,6]))) { error(); }
  414. if(any (any (2-a != [1,0;-1,-2]))) { error(); }
  415. if(any (any (a-2 != [-1,0;1,2]))) { error(); }
  416. if(any (any (2*a != [2,4;6,8]))) { error(); }
  417. if(any (any ((a./2 != [0.5,1;1.5,2])))) { error(); }
  418. if(any (any (a*a != [7,10;15,22]))) { error(); }
  419. if(any (any (a*a*a != [37,54;81,118]))) { error(); }
  420. if(any (any (a .* a != [1,4;9,16]))) { error(); }
  421. if(any (any (a./a != [1,1;1,1]))) { error(); }
  422. if(any (any (a' != [1,3;2,4]))) { error(); }
  423.  
  424. if(any(any(rand(3,3)^0 != eye(3,3)))) { error(); }
  425. if(any(any(rand(3,3).^0 != ones(3,3)))) { error(); }
  426. if(any(any(rand(1,3).^0 != ones(1,3)))) { error(); }
  427. if(any(any(rand(3,1).^0 != ones(3,1)))) { error(); }
  428. if(any(any(1.^zeros(3,1) != ones(3,1)))) { error(); }
  429. if(any(any(1.^zeros(1,3) != ones(1,3)))) { error(); }
  430.  
  431. if(any ([1;2;3]+[4;5;6] != [5;7;9])) {
  432.   error();
  433. }
  434. if(any ([1;2;3]-[4;5;6] != [-3;-3;-3])) {
  435.   error();
  436. }
  437. if(any ([2;2;2] ./ [1;1;1] != [2;2;2])) {
  438.   error();
  439. }
  440. if(any ([1;2;3] .* [4;5;6] != [4;10;18])) {
  441.   error();
  442. }
  443.  
  444. //
  445. // Test row-wise matrix addition
  446. //
  447.  
  448. a = [1,2,3]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  449. if (!all (all (a .+ b == [2,4,6;5,7,9;8,10,12;11,13,15]))) { error (); }
  450. if (!all (all (b .+ a == [2,4,6;5,7,9;8,10,12;11,13,15]))) { error (); }
  451.  
  452. a = [1,2,3] + [1,2,3]*1i; 
  453. b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
  454. c = [2,4,6;5,7,9;8,10,12;11,13,15] + [2,4,6;5,7,9;8,10,12;11,13,15]*1i;
  455. if (!all (all (a .+ b == c))) { error (); }
  456. if (!all (all (b .+ a == c))) { error (); }
  457.  
  458. printf("\t\tpassed matrix row-wise add test...\n");
  459.  
  460. //
  461. // Test row-wise matrix subtraction
  462. //
  463.  
  464. a = [1,1,1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  465. if (!all (all (a .- b == -[0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  466. if (!all (all (b .- a ==  [0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  467.  
  468. a = [1,1,1] + [1,1,1]*1i;
  469. b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
  470. c = [0,1,2;3,4,5;6,7,8;9,10,11] + [0,1,2;3,4,5;6,7,8;9,10,11]*1i;
  471. if (!all (all (a .- b == -c))) { error (); }
  472. if (!all (all (b .- a ==  c))) { error (); }
  473.  
  474. printf("\t\tpassed matrix row-wise subtraction test...\n");
  475.  
  476. //
  477. // Test col-wise matrix addition
  478. //
  479.  
  480. a = [1;1;1;1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  481. if (!all (all (a .+ b == [2,3,4;5,6,7;8,9,10;11,12,13]))) { error (); }
  482. if (!all (all (b .+ a == [2,3,4;5,6,7;8,9,10;11,12,13]))) { error (); }
  483.  
  484. a = [1;1;1;1] + [1;1;1;1]*1i;
  485. b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
  486. c = [2,3,4;5,6,7;8,9,10;11,12,13] + [2,3,4;5,6,7;8,9,10;11,12,13]*1i;
  487. if (!all (all (a .+ b == c))) { error (); }
  488. if (!all (all (b .+ a == c))) { error (); }
  489.  
  490. printf("\t\tpassed matrix col-wise add test...\n");
  491.  
  492. //
  493. // Test col-wise matrix subtraction
  494. //
  495.  
  496. a = [1;1;1;1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  497. if (!all (all (a .- b == -[0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  498. if (!all (all (b .- a ==  [0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  499.  
  500. a = [1;1;1;1] + [1;1;1;1]*1i;
  501. b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
  502. c = [0,1,2;3,4,5;6,7,8;9,10,11] + [0,1,2;3,4,5;6,7,8;9,10,11]*1i;
  503. if (!all (all (a .- b == -c))) { error (); }
  504. if (!all (all (b .- a ==  c))) { error (); }
  505.  
  506. printf("\t\tpassed matrix col-wise subtraction test...\n");
  507.  
  508. a = [1,2,3];
  509. b = [1,2,3;4,5,6;7,8,9];
  510. c = [1,4,9;4,10,18;7,16,27];
  511.  
  512. if (!all (all (a .* b == c))) { error (); }
  513. if (!all (all (b .* a == c))) { error (); }
  514.  
  515. za = a + rand (size (a))*1j;
  516. zb = b + rand (size (b))*1j;
  517.  
  518. if (!all (all (za .* zb == [za;za;za] .* zb))) { error (); }
  519. if (!all (all (zb .* za == zb .* [za;za;za]))) { error (); }
  520. if (!all (all (a .* zb == [a;a;a] .* zb))) { error (); }
  521. if (!all (all (zb .* a == zb .* [a;a;a]))) { error (); }
  522. if (!all (all (za .* b == [za;za;za] .* b))) { error (); }
  523. if (!all (all (b .* za == b .* [za;za;za]))) { error (); }
  524.  
  525. printf("\t\tpassed matrix row-wise multiplication test...\n");
  526.  
  527. a = [1,2,3];
  528. b = [1,2,3;4,6,6;7,8,9];
  529. c = [1,1,1;4,3,2;7,4,3];
  530.  
  531. if (!all (all (b ./ a == c))) { error (); }
  532. if (!all (all ([a;a;a] ./ b == a ./ b))) { error (); }
  533. if (!all (all (b ./ [a;a;a] == b ./ a))) { error (); }
  534.  
  535. za = a + rand (size (a))*1j;
  536. zb = b + rand (size (b))*1j;
  537.  
  538. if (!all (all ([za;za;za] ./ zb == za ./ zb))) { error (); }
  539. if (!all (all (zb ./ [za;za;za] == zb ./ za))) { error (); }
  540. if (!all (all ([a;a;a] ./ zb == a ./ zb))) { error (); }
  541. if (!all (all (zb ./ [a;a;a] == zb ./ a))) { error (); }
  542. if (!all (all ([za;za;za] ./ b == za ./ b))) { error (); }
  543. if (!all (all (b ./ [za;za;za] == b ./ za))) { error (); }
  544.  
  545. printf("\t\tpassed matrix row-wise division test...\n");
  546.  
  547. a = [1;2;3];
  548. b = [1,2,3;4,5,6;7,8,9];
  549.  
  550. if (!all (all (a .* b == [a,a,a] .* b))) { error (); }
  551. if (!all (all (b .* a == b .* [a,a,a]))) { error (); }
  552.  
  553. za = a + rand (size (a))*1j;
  554. zb = b + rand (size (b))*1j;
  555.  
  556. if (!all (all (za .* zb == [za,za,za] .* zb))) { error (); }
  557. if (!all (all (zb .* za == zb .* [za,za,za]))) { error (); }
  558. if (!all (all (za .* b == [za,za,za] .* b))) { error (); }
  559. if (!all (all (b .* za == b .* [za,za,za]))) { error (); }
  560. if (!all (all (a .* zb == [a,a,a] .* zb))) { error (); }
  561. if (!all (all (zb .* a == zb .* [a,a,a]))) { error (); }
  562.  
  563. printf("\t\tpassed matrix column-wise multiplication test...\n");
  564.  
  565. a = [1;2;3];
  566. b = [1,2,3;4,6,6;7,8,9];
  567.  
  568. if (!all (all ([a,a,a] ./ b == a ./ b))) { error (); }
  569. if (!all (all (b ./ [a,a,a] == b ./ a))) { error (); }
  570.  
  571. za = a + rand (size(a))*1j;
  572. zb = b + rand (size(b))*1j;
  573.  
  574. if (!all (all ([za,za,za] ./ zb == za ./ zb))) { error (); }
  575. if (!all (all (zb ./ [za,za,za] == zb ./ za))) { error (); }
  576. if (!all (all ([za,za,za] ./ b == za ./ b))) { error (); }
  577. if (!all (all (b ./ [za,za,za] == b ./ za))) { error (); }
  578. if (!all (all ([a,a,a] ./ zb == a ./ zb))) { error (); }
  579. if (!all (all (zb ./ [a,a,a] == zb ./ a))) { error (); }
  580.  
  581. printf("\t\tpassed matrix column-wise division test...\n");
  582.  
  583.  
  584. //
  585. //--------------------- Test COMPLEX MATRICES --------------------------
  586. //
  587. //  Automatic creation, extension
  588. //
  589. printf("\t\tcomplex-matrices\n");
  590. if(any (any ((mm[3;3]=10+10i) != [0,0,0;0,0,0;0,0,10+10i]))) { error(); }
  591. a=[1,2,3;4,5,6];
  592. if(any (any ((a[3;1]=10+10i) != [1,2,3;4,5,6;10+10i,0,0]))) { error(); }
  593. //
  594. a = m3;
  595. if(any (any (a+a != [2+4i,4+6i;6+8i,10+12i]))) { error(); }
  596. if(any (any (a-a != zeros(2,2)))) { error(); }
  597. if(any (any (2+a != [3+2i,4+3i;5+4i,7+6i]))) { error(); }
  598. if(any (any (2-a != [1-2i,0-3i;-1-4i,-3-6i]))) { error(); }
  599. if(any (any (a-2 != [-1+2i,0+3i;1+4i,3+6i]))) { error(); }
  600. if(any (any (2*a != [2+4i,4+6i;6+8i,10+12i]))) { error(); }
  601. if(any (any (a./2 != [.5+1i,1+1.5i;1.5+2i,2.5+3i]))) { error(); }
  602. if(any (any (a*a != [-9+21i,-12+34i;-14+48i,-17+77i]))) { error(); }
  603. if(any (any (a*a*a != [-223+57i,-345+113i;-469+183i,-719+337i]))) { error(); }
  604. if(any (any (a .* a != [-3+4i,-5+12i;-7+24i,-11+60i]))) { error(); }
  605. //
  606. // The following test may not work on some machines
  607. //
  608. if(any (any (a./a != [1,1;1,1]))) { 
  609.   printf("\t\t***complex division inaccuracy, check manually***\n");
  610. }
  611.  
  612. if(any (any (a' != conj([1+2i,3+4i;2+3i,5+6i])))) { error(); }
  613. //  
  614. //--------------------- Test NULL MATRICES -------------------------
  615. //
  616. printf("\t\tnull-matrices\n");
  617. // Create a NULL matrix
  618. mnull = matrix();
  619. // Test it with SCALARS
  620. if( any([1,mnull] != 1)) {
  621.   error();
  622. }
  623. if( any([mnull,1] != 1)) {
  624.   error();
  625. }
  626. // Test with MATRIX construction
  627. m = [1,2;3,4;5,6];
  628. if( any([mnull;1] != [1])) {
  629.   error();
  630. }
  631. if( any([1;mnull] != [1])) {
  632.   error();
  633. }
  634. if( any([mnull;1,2,3] != [1,2,3])) {
  635.   error();
  636. }
  637. if( any([1,2,3;mnull] != [1,2,3])) {
  638.   error();
  639. }
  640. if(any (any ([mnull,m] != m))) {
  641.   error();
  642. }
  643. if(any (any ([m,mnull] != m))) {
  644.   error();
  645. }
  646. if(any (any ([mnull;m] != m))) {
  647.   error();
  648. }
  649. if(any (any ([m;mnull] != m))) {
  650.   error();
  651.  
  652. mnull = matrix();
  653. mnull[1] = [1];
  654. }
  655.  
  656. //--------------------- Test Matrix Multiply --------------------------
  657.  
  658. i = sqrt(-1);
  659. a = [1,2,3;4,5,6;7,8,9];
  660. b = [4,5,6;7,8,9;10,11,12];
  661. c = [ 48,  54,  60 ;
  662.      111, 126, 141 ;
  663.      174, 198, 222 ];
  664.  
  665. if (any (any (c != a*b))) { error ("failed Real-Real Multiply"); }
  666.  
  667. az = a + b*i;
  668. bz = b + a*i;
  669.  
  670. cz = [-18+141*i , -27+162*i , -36+183*i ;
  671.         9+240*i ,   0+279*i ,  -9+318*i ;
  672.        36+339*i ,  27+396*i ,  18+453*i ];
  673.  
  674. czz = [ 48+30*i ,  54+36*i  ,  60+42*i ;
  675.        111+66*i , 126+81*i  , 141+96*i ;
  676.        174+102*i, 198+126*i , 222+150*i ];
  677.  
  678. czzz = [ 48+111*i ,  54+126*i ,  60+141*i ;
  679.         111+174*i , 126+198*i , 141+222*i ;
  680.         174+237*i , 198+270*i , 222+303*i ];
  681.  
  682. if (any (any (cz != az*bz)))  { error ("failed Complex-Complex Multiply"); }
  683. if (any (any (czz != a*bz)))  { error ("failed Real-Complex Multiply"); }
  684. if (any (any (czzz != az*b))) { error ("failed Complex-Real Multiply"); }
  685.  
  686. a = [a,a];
  687. b = [b;b];
  688. c = [  96 , 108 , 120 ;
  689.       222 , 252 , 282 ;
  690.       348 , 396 , 444 ];
  691.  
  692. if (any (any (c != a*b))) { error ("failed Real-Real Multiply"); }
  693.  
  694. az = [az,az];
  695. bz = [bz;bz];
  696.  
  697. cz = [  -36+282*i ,  -54+324*i ,  -72+366*i ;
  698.          18+480*i ,    0+558*i ,  -18+636*i ;
  699.          72+678*i ,   54+792*i ,   36+906*i ];
  700.  
  701. czz = [  96+60*i  , 108+72*i  , 120+84*i  ;
  702.         222+132*i , 252+162*i , 282+192*i ;
  703.         348+204*i , 396+252*i , 444+300*i ];
  704.  
  705. czzz = [  96+222*i , 108+252*i , 120+282*i ;
  706.          222+348*i , 252+396*i , 282+444*i ;
  707.          348+474*i , 396+540*i , 444+606*i ];
  708.  
  709. if (any (any (cz != az*bz)))  { error ("failed Complex-Complex Multiply"); }
  710. if (any (any (czz != a*bz)))  { error ("failed Real-Complex Multiply"); }
  711. if (any (any (czzz != az*b))) { error ("failed Complex-Real Multiply"); }
  712.  
  713. printf("\t\tpassed matrix multiply test...\n");
  714.  
  715. //--------------------- Test STRING MATRICES --------------------------
  716. //
  717. printf("\t\tstring-matrices\n");
  718. sm = ["s1","sm2","sm3";"x1","x2","xxx3";"y1","y2","yyy3"];
  719. if(sm[1] != "s1") { error(); }
  720. if( sm[1;3] != "sm3" ) { error(); }
  721. if(any(sm[2,3;3] != ["xxx3";"yyy3"]) ) { error(); }
  722. if(any (any ((sm[1;1]="xx")!=["xx","sm2","sm3";"x1","x2","xxx3";"y1","y2","yyy3"])))
  723. {
  724.   error();
  725. }
  726. if( ((strm[1] = "strm")[1]) != "strm" ) { error(); }
  727.  
  728. //  Test string-matrix add functionality
  729.  
  730. sm = sm[1,2;1,2];
  731. if (any (any (("1_"+sm+"_2") != ["1_xx_2","1_sm2_2";"1_x1_2","1_x2_2"]))) {error();}
  732.  
  733. if ("c"+["1"] != "c1") { error (); }
  734. if (["c"]+"1" != "c1") { error (); }
  735.  
  736. //
  737. //---------------------------- List Tests --------------------------
  738. //
  739. //  List creation
  740. listest = << << 11; 12 >>; << 21; 22>> >>;
  741. if( listest.[1].[2] != 12 ) { error(); }
  742. if(any(<<a=10;b=1:4;c=[1,2,3;4,5,6;7,8,9]>>.b != [1,2,3,4])) { error(); }
  743. mlist.[0] = m;
  744. if(any(any(mlist.[0] != m))) { error(); }
  745. printf("\tlist tests\n");
  746.  
  747. //
  748. // Reset random number generator seed...
  749. //
  750.  
  751. rand("default");
  752.  
  753. //
  754. //-------------------------- Test abs () ------------------------------
  755. //
  756.  
  757. A = rand(5,5);
  758. T = ( A == abs (A));
  759. if (!all (all (A))) { error ("abs() incorrect"); }
  760. printf("\tabs()");
  761.  
  762. //
  763. //-------------------------- Test max () ------------------------------
  764. //
  765.  
  766. A = [1,10,100;2,20,200;1,2,3];
  767. B = A./2;
  768. ZA = A + rand (3,3)*A*1i;
  769. ZB = B + rand (3,3)*B*1i;
  770. if (!all (max (A) == [2,20,200])) { error( "max() incorrect"); }
  771. if (max (max(A)) != 200) { error ("max() incorrect"); }
  772. if (any (any (max (A, B) != A))) { error (); }
  773. if (any (any (max (B, A) != max (A, B)))) { error (); }
  774. if (any (any (max (ZB, ZA) != max (ZA, ZB)))) { error (); }
  775. if (any (any (max (B, ZA) != max (ZA, B)))) { error (); }
  776. if (any (any (max (ZB, A) != max (A, ZB)))) { error (); }
  777. printf("\tmax()");
  778.  
  779. //
  780. //-------------------------- Test min () ------------------------------
  781. //
  782.  
  783. if (!all (min (A) == [1,2,3])) { error( "min() incorrect"); }
  784. if (min (min(A)) != 1) { error ("min() incorrect"); }
  785. if (any (any (min (A, B) != B))) { error (); }
  786. if (any (any (min (B, A) != min (A, B)))) { error (); }
  787. if (any (any (min (ZB, ZA) != min (ZA, ZB)))) { error (); }
  788. if (any (any (min (B, ZA) != min (ZA, B)))) { error (); }
  789. if (any (any (min (ZB, A) != min (A, ZB)))) { error (); }
  790. printf("\tmin()");
  791.  
  792. //
  793. //-------------------------- Test maxi () -----------------------------
  794. //
  795.  
  796. if (!all (maxi (A) == [2,2,2])) { error( "maxi() incorrect"); }
  797. if (maxi (maxi(A)) != 1) { error ("maxi() incorrect"); }
  798. printf("\tmaxi()");
  799.  
  800. //
  801. //-------------------------- Test mini () -----------------------------
  802. //
  803.  
  804. if (!all (mini (A) == [1,3,3])) { error( "mini() incorrect"); }
  805. if (mini (mini(A)) != 1) { error ("mini() incorrect"); }
  806. printf("\tmini()");
  807.  
  808. //
  809. //-------------------------- Test floor () ----------------------------
  810. //
  811.  
  812. if (floor (1.9999) != 1) { error ("floor() output incorrect"); }
  813. if (!all (all (floor ([1.99,1.99;2.99,2.99]) == [1,1;2,2]))) {
  814.   error ("floor output incorrect");
  815. }
  816. printf("\tfloor()");
  817.  
  818. //
  819. //-------------------------- Test ceil () 0----------------------------
  820. //
  821.  
  822. if (ceil (1.9999) != 2) { error ("ceil() output incorrect"); }
  823. if (!all (all (ceil ([1.99,1.99;2.99,2.99]) == [2,2;3,3]))) {
  824.   error ("ceil output incorrect");
  825. }
  826. printf("\tceil()\n");
  827.  
  828. //
  829. //-------------------------- Test round () ----------------------------
  830. //
  831.  
  832. if (round (1.8) != 2) { error ("round() output incorrect"); }
  833. if (round (1.4) != 1) { error ("round() output incorrect"); }
  834. if (!all (all (round ([1.99,1.99;2.4,2.4]) == [2,2;2,2]))) {
  835.   error ("round output incorrect");
  836. }
  837. printf("\tround()");
  838.  
  839. //
  840. //-------------------------- Test sum () ------------------------------
  841. //
  842.  
  843. S = [1:4; 4:7; 8:11];
  844. if (sum (S[1;]) != 10) { error ("sum() incorrect"); }
  845. if (sum (S[3;]) != 38) { error ("sum() incorrect"); }
  846. if (!all (all (sum (S) == [13,16,19,22]))) { error ("sum() incorrect"); }
  847. printf("\tsum()");
  848.  
  849. //
  850. //-------------------------- Test int () ------------------------------
  851. //
  852.  
  853. if (int (1.9999) != 1) { error ("int() output incorrect"); }
  854. if (!all (all (int ([1.99,1.99;2.99,2.99]) == [1,1;2,2]))) {
  855.   error ("int() output incorrect");
  856. }
  857. printf("\tint()");
  858.  
  859. //
  860. //-------------------------- Test mod () ------------------------------
  861. //
  862.  
  863. if (mod (1,1) != 0) { error ("mod() output incorrect"); }
  864. if (mod (4,2) != 0) { error ("mod() output incorrect"); }
  865. if (mod (3,2) != 1) { error ("mod() output incorrect"); }
  866. if (mod (5,3) != 2) { error ("mod() output incorrect"); }
  867. printf("\tmod()");
  868.  
  869. //
  870. //-------------------------- Test find () ------------------------------
  871. //
  872.  
  873. if (find ([0,1]) != 2) { error ("find() output incorrect"); }
  874. if (find ([1,0]) != 1) { error ("find() output incorrect"); }
  875. if (find ([0,1+1i]) != 2) { error ("find() output incorrect"); }
  876. if (find ([1+0i,0]) != 1) { error ("find() output incorrect"); }
  877. if (find ([0+1i,0]) != 1) { error ("find() output incorrect"); }
  878.  
  879. printf("\tfind()\n");
  880.  
  881. //
  882. //-------------------------- Test rand () -----------------------------
  883. //
  884. rand ("normal", 5, 1);
  885. xrand = rand(4000,1);
  886.  
  887. mean = function(x)
  888. {
  889.   local(m);
  890.  
  891.   m = size (x)[1];
  892.   if( m == 1 ) 
  893.   { 
  894.     m = size (x)[2];
  895.   }
  896.  
  897.   return sum( x ) / m;
  898. };
  899.  
  900. std = function(x)
  901. {
  902.   local(i, m, s);
  903.  
  904.   if(class(x) != "num") { error("std() requires NUMERICAL input"); }
  905.  
  906.   m = x.nr;
  907.   if( m == 1 ) 
  908.   { 
  909.     return sqrt( sum( (x - mean(x)) .^ 2 ) / (x.nc - 1) );
  910.   else
  911.     for( i in 1:x.nc) {
  912.       s[i] = sqrt( sum( (x[;i] - mean(x[;i])) .^ 2 ) / (x.nr - 1) );
  913.     }
  914.     return s;
  915.   }
  916. };
  917.  
  918. if (!(mean (xrand) > 4.9 && mean (xrand) < 5.1)) 
  919.   { error ("error in random"); }
  920. if (!(std (xrand) > 0.9 && std (xrand) < 1.1))
  921.   { error ("error in random"); }
  922. printf("\tpassed rand test...\n");
  923.  
  924. rand("default");
  925.  
  926. //
  927. //-------------------------- Test norm () -----------------------------
  928. //
  929.  
  930. tn = [1,2,3,4;2,1,2,3;3,2,1,2;4,3,2,1  ];
  931. if (norm(tn,"m") != 4 ) { error ("incorrect norm computation"); }
  932. if (norm(tn,"1") != 10) { error ("incorrect norm computation"); }
  933. if (norm(tn,"i") != 10) { error ("incorrect norm computation"); }
  934. printf("\tpassed norm test...\n");
  935.  
  936. //
  937. //-------------------------- Test qr () ------------------------------
  938. //
  939.  
  940. a = rand (4, 4);
  941. qa = qr (a);
  942. if (max (max (abs (qa.q*qa.r - a)))/(X*norm (a)*a.nr) > eps)
  943.   { error ("possible qr() problems"); }
  944.  
  945. z = rand (4,4) + rand(4,4)*1i;
  946. qz = qr (z);
  947. if (max (max (abs (qz.q*qz.r -  z)))/(X*norm (z)*z.nr) > eps)
  948.   { error ("possible qr() problems"); }
  949.  
  950. printf("\tpassed qr test...\n");
  951.  
  952. //
  953. //-------------------------- Test schur () ----------------------------
  954. //
  955.  
  956. a = rand (4, 4);
  957. sa = schur (a);
  958. if (max (max (abs (sa.z*sa.t*sa.z' - a)))/(X*norm (a)*a.nr) > eps)
  959.   { error ("possible schur() problems"); }
  960.  
  961. z = rand (4,4) + rand(4,4)*1i;
  962. sz = schur (z);
  963. if (max (max (abs (sz.z*sz.t*sz.z' - z)))/(X*norm (z)*z.nr) > eps)
  964.   { error ("possible schur() problems"); }
  965.  
  966. printf("\tpassed schur test...\n");
  967.  
  968. //
  969. //-------------------------- Test schord () ----------------------------
  970. //
  971.  
  972. s2a = schord (sa, 2, 4);
  973. if (max (max (abs (s2a.z*s2a.t*s2a.z' - a)))/(X*norm (a)*a.nr) > eps)
  974.   { error ("possible schord() problems"); }
  975.  
  976. s2z = schord (sz, 3, 1);
  977. if (max (max (abs (s2z.z*s2z.t*s2z.z' - z)))/(X*norm (z)*z.nr) > eps)
  978.   { error ("possible schord() problems"); }
  979.  
  980. printf("\tpassed schord test...\n");
  981.  
  982.  
  983. //
  984. //-------------------------- Test chol () -----------------------------
  985. //
  986.  
  987. c = rand(4,4);
  988. c = symm (c + eye (size (c))*norm(c)*c.nr);
  989. u = chol (c);
  990. if (max (max (abs (u'*u - c)))/(X*norm (c)*c.nr) > eps)
  991.   { error ("possible chol() problems"); }
  992.  
  993. cz = rand(4,4);
  994. cz = symm (cz + eye (size (cz))*X*norm(cz)*cz.nr);
  995. uz = chol (cz);
  996. if (max (max (abs (uz'*uz - cz)))/(X*norm (cz)*cz.nr) > eps)
  997.   { error ("possible chol() problems"); }
  998.  
  999. printf("\tpassed chol test...\n");
  1000.  
  1001. //
  1002. //-------------------------- Test solve () -----------------------------
  1003. //
  1004.  
  1005. a = rand(4,4);
  1006. while (1/rcond (a) > 100) {  a = rand(4,4); }
  1007. b = ones(4,1);
  1008. x = solve (a,b);
  1009. if (max (max (abs(a*x - b)))/(X*norm (a)*a.nr) > eps)
  1010.   { 
  1011.     printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  1012.     printf ("\tA*X - B:\n");
  1013.     abs (a*x - b)
  1014.     error ("possible solve() problems\n");
  1015.   }
  1016.  
  1017. az = rand(4,4) + rand(4,4)*1j;
  1018. while (1/rcond (az) > 100) { az = rand(4,4) + rand(4,4)*1j; }
  1019. bz = rand(4,1) + rand(4,1)*1j;
  1020. xz = solve (az,bz);
  1021. if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
  1022.   { 
  1023.     printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  1024.     printf ("\tA*X - B:\n");
  1025.     abs (az*xz - bz)
  1026.     error ("possible solve() problems\n");
  1027.   }
  1028.  
  1029. printf("\tpassed solve test...\n");
  1030.  
  1031.  
  1032. //
  1033. //-------------------------- Test lu() / factor() -----------------------------
  1034. //
  1035.  
  1036. tril = function(x, k) 
  1037. {
  1038.   local(i, j, nr, nc, y);
  1039.  
  1040.   if (!exist (k)) { k = 0; }
  1041.   nr = x.nr; nc = x.nc;
  1042.   if(k > 0) 
  1043.   { 
  1044.     if (k > (nc - 1)) { error ("tril: invalid value for k"); }
  1045.   else
  1046.     if (abs (k) > (nr - 1)) { error ("tril: invalid value for k"); }
  1047.   }
  1048.  
  1049.   y = zeros(nr, nc);
  1050.  
  1051.   for(i in max( [1,1-k] ):nr) {
  1052.     j = 1:min( [nc, i+k] );
  1053.     y[i;j] = x[i;j];
  1054.   }
  1055.  
  1056.   return y;
  1057. };
  1058.  
  1059. triu = function(x, k) 
  1060. {
  1061.   local(i, j, nr, nc, y);
  1062.  
  1063.   if (!exist (k)) { k = 0; }
  1064.   nr = x.nr; nc = x.nc;
  1065.  
  1066.   if(k > 0) 
  1067.   { 
  1068.     if (k > (nc - 1)) { error ("triu: invalid value for k"); }
  1069.   else
  1070.     if (abs (k) > (nr - 1)) { error ("triu: invalid value for k"); }
  1071.   }
  1072.  
  1073.   y = zeros(nr, nc);
  1074.  
  1075.   for(j in max( [1,1+k] ):nc) {
  1076.     i = 1:min( [nr, j-k] );
  1077.     y[i;j] = x[i;j];
  1078.   }
  1079.  
  1080.   return y;
  1081. };
  1082.  
  1083. static (swap);
  1084.  
  1085. lu = function ( A )
  1086. {
  1087.   local (i, l, u, pvt, x)
  1088.  
  1089.   if (A.nr != A.nc) { error ("lu() requires square A"); }
  1090.  
  1091.   x = factor (A);    // Do the factorization
  1092.  
  1093.   //
  1094.   // Now create l, u, and pvt from lu and pvt.
  1095.   //
  1096.  
  1097.   l = tril (x.lu, -1) + eye (size (x.lu));
  1098.   u = triu (x.lu);
  1099.   pvt = eye (size (x.lu));
  1100.  
  1101.   //
  1102.   // Now re-arange the columns of pvt
  1103.   //
  1104.  
  1105.   for (i in 1:max (size (x.lu)))
  1106.   {
  1107.     pvt = pvt[ ; swap (1:pvt.nc, i, x.pvt[i]) ];
  1108.   }
  1109.   return << l = l; u = u; pvt = pvt >>;
  1110. };
  1111.  
  1112. //
  1113. //  In vector V, swap elements I, J
  1114. //
  1115.  
  1116. swap = function ( V, I, J )
  1117. {
  1118.   local (v, tmp);
  1119.   v = V;
  1120.   tmp = v[I];
  1121.   v[I] = v[J];
  1122.   v[J] = tmp;
  1123.   return v;
  1124. };
  1125.  
  1126.  
  1127. a = rand(4,4);
  1128. while (1/rcond (a) > 100) {  a = rand(4,4); }
  1129. lua = lu (a);
  1130. if (max (max (abs(a - lua.pvt*lua.l*lua.u)))/(X*norm (a)*a.nr) > eps)
  1131.   { 
  1132.     printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  1133.     printf ("\tA - p*l*u:\n");
  1134.     abs (a - lua.pvt*lua.l*lua.u)
  1135.     error ("possible lu()/factor() problems\n");
  1136.   }
  1137.  
  1138. az = rand(4,4) + rand(4,4)*1j;
  1139. while (1/rcond (az) > 100) { az = rand(4,4) + rand(4,4)*1j; }
  1140. luz = lu (az);
  1141. if (max (max (abs (az - luz.pvt*luz.l*luz.u)))/(X*norm (az)*az.nr) > eps)
  1142.   { 
  1143.     printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  1144.     printf ("\tA - p*l*u:\n");
  1145.     abs (az - luz.ovt*luz.l*luz.u)
  1146.     error ("possible lu()/factor()() problems\n");
  1147.   }
  1148.  
  1149. printf("\tpassed lu/factor test...\n");
  1150.  
  1151.  
  1152. //
  1153. //-------------------------- Test svd ()   -----------------------------
  1154. //
  1155.  
  1156. a = rand(4,4);
  1157. s = svd (a);
  1158. if (max (max (abs (s.u*diag(s.sigma)*s.vt - a)))/(X*norm (a)*a.nr) > eps)
  1159.   { error ("possible svd() problems"); }
  1160.  
  1161. z = rand(4,4) + rand(4,4)*1j;
  1162. sz = svd (z);
  1163. if (max (max (abs (sz.u*diag(sz.sigma)*sz.vt - z)))/(X*norm (z)*z.nr) > eps)
  1164.   { error ("possible svd() problems"); }
  1165.  
  1166. printf("\tpassed svd test...\n");
  1167.  
  1168.  
  1169. //
  1170. //-------------------------- Test hess ()  -----------------------------
  1171. //
  1172.  
  1173. a = rand(4,4);
  1174. h = hess (a);
  1175. if (max (max (abs (h.p*h.h*h.p' - a)))/(X*norm (a)*a.nr) > eps)
  1176.   { error ("possible hess() problems"); }
  1177.  
  1178. z = rand(4,4) + rand(4,4)*1j;
  1179. hz = hess (z);
  1180. if (max (max (abs (hz.p*hz.h*hz.p' - z)))/(X*norm (z)*z.nr) > eps)
  1181.   { error ("possible hess() problems"); }
  1182.  
  1183. printf("\tpassed hess test...\n");
  1184.  
  1185. //
  1186. //-------------------------- Test lyap () ------------------------------
  1187. //
  1188.  
  1189. lyap = function ( A, B, C )
  1190. {
  1191.   local (X, sa, sb, tc)
  1192.  
  1193.   if (!exist (B)) 
  1194.   { 
  1195.     B = A';    // Solve the special form: A*X + X*A' = -C
  1196.   }
  1197.  
  1198.   if ((A.nr != A.nc) || (B.nr != B.nc) || (C.nr != A.nr) || (C.nc != B.nr)) {
  1199.     error ("Dimensions do not agree.");
  1200.   }
  1201.  
  1202.   //
  1203.   // Schur decomposition on A and B
  1204.   //
  1205.  
  1206.   sa = schur (A);
  1207.   sb = schur (B);
  1208.  
  1209.   //
  1210.   // transform C
  1211.   //
  1212.  
  1213.   tc = sa.z' * C * sb.z;
  1214.  
  1215.   X = sylv (sa.t, sb.t, tc);
  1216.  
  1217.   //
  1218.   // Undo the transformation
  1219.   //
  1220.  
  1221.   X = sa.z * X * sb.z';
  1222.  
  1223.   return X;
  1224. };
  1225.  
  1226. a = rand (5,5);
  1227. b = rand (5,5);
  1228. c = rand (5,5);
  1229. while (1/rcond (a) > 100) {  a = rand(5,5); }
  1230. while (1/rcond (b) > 100) {  b = rand(5,5); }
  1231. while (1/rcond (c) > 100) {  c = rand(5,5); }
  1232.  
  1233. x = lyap (a, b, c);
  1234. if (max (max (abs (a*x + x*b + c)))/(X*norm(c)*norm(a)*norm(b)) > eps)
  1235.   { error ("possible problems with lyap() or sylv()"); }
  1236.  
  1237. printf("\tpassed lyap test...\n");
  1238.  
  1239. //
  1240. //-------------------------- Test eig () ------------------------------
  1241. //
  1242.  
  1243. trace = function(m) 
  1244. {
  1245.   local(i, tr);
  1246.  
  1247.   if(m.class != "num") { 
  1248.     error("must provide NUMERICAL input to trace()");
  1249.   }
  1250.  
  1251.   tr = 0;
  1252.   for(i in 1:min( [m.nr, m.nc] )) {
  1253.     tr = tr + m[i;i];
  1254.   }
  1255.  
  1256.   return tr;
  1257. };
  1258.  
  1259. eye = function( m , n ) 
  1260. {
  1261.   local(i, N, new);
  1262.  
  1263.   if (!exist (n))
  1264.   {
  1265.     if(m.n != 2) { error("only 2-el MATRIX allowed as eye() arg"); }
  1266.     new = zeros (m[1], m[2]);
  1267.     N = min ([m[1], m[2]]);
  1268.   else
  1269.     if (class (m) == "string" || class (n) == "string") {
  1270.       error ("eye(), string arguments not allowed");
  1271.     }
  1272.     if (max (size (m)) == 1 && max (size (n)) == 1)
  1273.     {
  1274.       new = zeros (m[1], n[1]);
  1275.       N = min ([m[1], n[1]]);
  1276.     else
  1277.       error ("matrix arguments to eye() must be 1x1");
  1278.     }
  1279.   }
  1280.   for(i in 1:N)
  1281.   {
  1282.     new[i;i] = 1.0;
  1283.   }
  1284.   return new;
  1285. };
  1286.  
  1287. //
  1288. // Standard eigenvalue problem
  1289. //
  1290.  
  1291. a = rand(10,10);
  1292. sa = symm (a);
  1293. ta = trace (a);
  1294. tsa = trace (sa);
  1295.  
  1296. z = rand(10,10) + rand(10,10)*1i;
  1297. sz = symm (z);
  1298. tz = trace (z);
  1299. tsz = trace (sz);
  1300.  
  1301. tol = 1.e-7;
  1302.  
  1303. if (!(ta < sum(eig(a).val) + tol && ta > sum(eig(a).val) - tol))
  1304. {
  1305.   error ("error in eig");
  1306. }
  1307. if (!(tsa < sum(eig(sa).val) + tol && tsa > sum(eig(sa).val) - tol))
  1308. {
  1309.   error ("error in eig");
  1310. }
  1311. if (!(tz < sum(eig(z).val) + tol && tz > sum(eig(z).val) - tol))
  1312. {
  1313.   error ("error in eig");
  1314. }
  1315. if (!(tsz < sum(eig(sz).val) + tol && tsz > sum(eig(sz).val) - tol))
  1316. {
  1317.   error ("error in eig");
  1318. }
  1319.  
  1320. //
  1321. // Generalized eigenvalue problem
  1322. //
  1323.  
  1324. b = rand(10,10);
  1325. sb = symm (b) + eye(size(b))*3;
  1326. tb = trace (b);
  1327. tsb = trace (sb);
  1328.  
  1329. zb = rand(10,10) + rand(10,10)*1i;
  1330. szb = symm (zb) + eye(size(zb))*3;
  1331. tzb = trace (zb);
  1332. tszb = trace (szb);
  1333.  
  1334. eig(a,b);    // not sure of a good way to check these yet
  1335. eigs(sa,sb);
  1336. eigs(sa,sb);
  1337.  
  1338. eig(z, zb);
  1339. eigs(sz, szb);
  1340. eigs(sz, szb);
  1341.  
  1342. printf("\tpassed eig test...\n");
  1343.  
  1344. //
  1345. //-------------------------- Test fft () -----------------------------
  1346. //
  1347.  
  1348. if (100 != fft(ones(100,1))[1]) { error ("error in fft()"); }
  1349. printf("\tpassed fft test...\n");
  1350.  
  1351. //
  1352. //-------------------------- Test printf () --------------------------
  1353. //
  1354.  
  1355. sprintf (tmp, "%*.*d %*.*d %s %*.*f f\n", 5,3,2, 8,7,3, "string", 3, 4, 1234e-2);
  1356. if (!(tmp == "  002  0000003 string 12.3400 f\n")) { error ("sprintf() error"); }
  1357. printf("\tpassed sprintf test...\n");
  1358.  
  1359. //
  1360. //------------------------- Test strtod()  ----------------------------
  1361. //
  1362.  
  1363. if (123.456 != strtod ("123.456")) { error (); }
  1364. if (!all (all ([1,2;3,4] == strtod (["1","2";"3","4"]))))
  1365.   { error (); }
  1366. printf("\tpassed strtod test...\n");
  1367.  
  1368. //
  1369. //------------------------- Test getline()  ---------------------------
  1370. //
  1371. //
  1372.  
  1373. close( "test.getline" );
  1374.  
  1375. x = getline( "test.getline" );
  1376. if ( x.[1] !=  123.456 ) { error(); }
  1377. if ( x.[2] != -123.456 ) { error(); }
  1378. if ( x.[3] !=  123.456 ) { error(); }
  1379.  
  1380. x = getline( "test.getline" );
  1381. if ( x.[1] !=  .123 ) { error(); }
  1382. if ( x.[2] != -.123 ) { error(); }
  1383. if ( x.[3] !=  .123 ) { error(); }
  1384.  
  1385. x = getline( "test.getline" );
  1386. if ( x.[1] !=  123 ) { error(); }
  1387. if ( x.[2] != -123 ) { error(); }
  1388. if ( x.[3] !=  123 ) { error(); }
  1389.  
  1390. x = getline( "test.getline" );
  1391. if ( x.[1] !=  1e6 ) { error(); }
  1392. if ( x.[2] != -1e6 ) { error(); }
  1393. if ( x.[3] !=  1e6 ) { error(); }
  1394.  
  1395. x = getline( "test.getline" );
  1396. if ( x.[1] !=  1e5 ) { error(); }
  1397. if ( x.[2] != -1e5 ) { error(); }
  1398. if ( x.[3] !=  1e5 ) { error(); }
  1399.  
  1400. x = getline( "test.getline" );
  1401. if ( x.[1] !=  123.456e3 ) { error(); }
  1402. if ( x.[2] != -123.456e3 ) { error(); }
  1403. if ( x.[3] !=  123.456e3 ) { error(); }
  1404.  
  1405. x = getline( "test.getline" );
  1406. if ( x.[1] !=  123.456e3 ) { error(); }
  1407. if ( x.[2] != -123.456e3 ) { error(); }
  1408. if ( x.[3] !=  123.456e3 ) { error(); }
  1409.  
  1410. x = getline( "test.getline" );
  1411. if ( x.[1] !=  123.456e-3 ) { error(); }
  1412. if ( x.[2] != -123.456e-3 ) { error(); }
  1413. if ( x.[3] !=  123.456e-3 ) { error(); }
  1414.  
  1415. x = getline( "test.getline" );
  1416. if ( x.[1] !=  .123e3 ) { error(); }
  1417. if ( x.[2] != -.123e3 ) { error(); }
  1418. if ( x.[3] !=  .123e3 ) { error(); }
  1419.  
  1420. x = getline( "test.getline" );
  1421. if ( x.[1] !=  123e3 ) { error(); }
  1422. if ( x.[2] != -123e3 ) { error(); }
  1423. if ( x.[3] !=  123e3 ) { error(); }
  1424.  
  1425. x = getline( "test.getline" );
  1426. if ( x.[1] != "123abc" ) { error(); }
  1427. if ( x.[2] != "abc123" ) { error(); }
  1428. if ( x.[3] != "123.abc" ) { error(); }
  1429.  
  1430. x = getline( "test.getline" );
  1431. if ( x.[1] != "quoted string" ) { error(); }
  1432. if ( x.[2] != "q string with escapes \n \t \" " ) { error(); }
  1433.  
  1434. x = getline( "test.getline" );
  1435. if ( x.[1] != "quoted string" ) { error(); }
  1436. if ( x.[2] !=  1.23e3 ) { error(); }
  1437. if ( x.[3] !=  100 ) { error(); }
  1438. if ( x.[4] != "q string with escapes \n \t \" " ) { error(); }
  1439. if ( x.[5] !=  200 ) { error(); }
  1440.  
  1441. printf("\tpassed getline() test...\n");
  1442.  
  1443. //
  1444. //---------------------- Test readb()/writeb()  --------------------
  1445. //
  1446. a = rand (5,5);
  1447. z = rand(3,3) + rand(3,3)*1j;
  1448. strm = what()[1:5;1:5];
  1449. pi = 4*atan(1.0);
  1450. sc = 2*pi;
  1451. str = "this is a sample string\ttab";
  1452. l = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1453.  
  1454. writeb ("jnk.rb", a, z, strm, sc, str, l);
  1455.  
  1456. # Set aside matrices for later tests
  1457.  
  1458. check = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1459.  
  1460. clear (a, z, strm, sc, str, l);
  1461.  
  1462. close ("jnk.rb");
  1463.  
  1464. readb ("jnk.rb");
  1465.  
  1466. #
  1467. # Now do checks
  1468. #
  1469.  
  1470. if (!all (all (a == check.a))) { error (); }
  1471. if (!all (all (z == check.z))) { error (); }
  1472. if (!all (all (strm == check.strm))) { error (); }
  1473. if (sc != check.sc) { error (); }
  1474. if (str != check.str) { error (); }
  1475.  
  1476. if (length (l) != 5) { error (); }
  1477. if (!all (all (l.a == check.a))) { error (); }
  1478. if (!all (all (l.z == check.z))) { error (); }
  1479. if (!all (all (l.strm == check.strm))) { error (); }
  1480. if (l.sc != check.sc) { error (); }
  1481. if (l.str != check.str) { error (); }
  1482.  
  1483.  
  1484. printf("\tpassed binary I/O test...\n");
  1485.  
  1486. //
  1487. //---------------------- Test read()/write()  --------------------
  1488. //
  1489. a = rand (5,5);
  1490. z = rand(3,3) + rand(3,3)*1j;
  1491. strm = what()[1:5;1:5];
  1492. pi = 4*atan(1.0);
  1493. sc = 2*pi;
  1494. str = "this is a sample string\ttab";
  1495. l = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1496.  
  1497. write ("jnk.ra", a, z, strm, sc, str, l);
  1498.  
  1499. # Set aside matrices for later tests
  1500.  
  1501. check = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1502.  
  1503. clear (a, z, strm, sc, str, l);
  1504.  
  1505. close ("jnk.ra");
  1506.  
  1507. read ("jnk.ra");
  1508.  
  1509. #
  1510. # Now do checks
  1511. #
  1512.  
  1513. if (a.nr != 5 || a.nc != 5) { error (); }
  1514. if (z.nr != 3 || z.nc != 3) { error (); }
  1515. if (strm.nr != 5 || strm.nc != 5) { error (); }
  1516. if (str != check.str) { error (); }
  1517.  
  1518. if (length (l) != 5) { error (); }
  1519. if (l.a.nr != 5 || l.a.nc != 5) { error (); }
  1520. if (l.z.nr != 3 || l.z.nc != 3) { error (); }
  1521. if (l.strm.nr != 5 || l.strm.nc != 5) { error (); }
  1522. if (l.str != check.str) { error (); }
  1523.  
  1524. printf("\tpassed ascii I/O test...\n");
  1525.  
  1526. //
  1527. //------------------------- Fibonacci Test -------------------------
  1528. //
  1529. //  Calculate Fibonacci numbers
  1530. //
  1531. i=1; 
  1532. while ( i < 2 ) { 
  1533.   i=i+1;
  1534.   a=0; b=1;
  1535.   while ( b < 10000 ) {
  1536.     c = b;
  1537.     b = a+b;
  1538.     a = c;
  1539.   }
  1540. }
  1541. if ( b != 10946 ) {
  1542.   error("failed fibonacci test");
  1543. else
  1544.   printf("\tpassed fibonacci test...\n");
  1545. }
  1546.  
  1547. //
  1548. //------------------------- Factorial Test -------------------------
  1549. //
  1550. fac = function(a) {
  1551.   if(a <= 1) {
  1552.     return 1
  1553.   else
  1554.     return a*fac(a-1)
  1555.   }
  1556. };
  1557. if(fac(10) != 3628800) { error(); else printf("\tpassed factorial test...\n"); }
  1558. //
  1559. //--------------------------- ACK Test ----------------------------
  1560. //
  1561. ack = function(a, b) {
  1562.   if(a == 0) { return b + 1; }
  1563.   if(b == 0) { return $self(a-1, 1); }
  1564.   return $self(a-1, $self(a, b-1));
  1565. };
  1566.  
  1567. if(ack(2,2) != 7) {
  1568.   error("error in ack() test");
  1569. else
  1570.   printf("\tpassed ACK test...\n");
  1571. }
  1572. //
  1573. //------------------------- Prime Test -----------------------------
  1574. //
  1575. // An example that finds all primes less than limit
  1576. //
  1577. primes = function (limit) {
  1578.   local(prime, cnt, i, j, k);
  1579.  
  1580.   i = 1; cnt = 0;
  1581.   for(k in 2:limit) {
  1582.     j = 2;
  1583.     while(mod(k,j) != 0) {
  1584.       j++;
  1585.     }
  1586.     if(j == k) {            // Found prime
  1587.       cnt++;
  1588.       prime[i;1] = k;
  1589.       i++;
  1590.     }
  1591.   }
  1592.   return prime;
  1593. };
  1594.  
  1595. if(max(size(primes(100))) != 25) { 
  1596.   error("error in prime test");
  1597. else
  1598.   printf("\tpassed prime test...\n");
  1599. }
  1600.  
  1601. //
  1602. //--------------------------- Fib Min Test -----------------------------
  1603. //
  1604. //    fibmin() will minimize an arbitrary function 
  1605. //    in 1D using Fibonacci search
  1606.  
  1607. f065 = function(x)
  1608. {
  1609.     return (x - 0.65) * (x - 0.65);
  1610. };
  1611.  
  1612. fib = function(x)
  1613. {
  1614.     local (n, a, b);
  1615.  
  1616.     a = 1;
  1617.     b = 1;
  1618.     if (x >= 2)
  1619.     {
  1620.         n = x - 1;
  1621.         for (n in n:1:-1)
  1622.         {
  1623.             c = b;
  1624.             b = a + b;
  1625.             a = c;
  1626.             n = n - 1;
  1627.         }
  1628.     }
  1629.     return b;
  1630. };
  1631.  
  1632. //  Minimize a 1D function using Fibonacci search
  1633. //  f = function to minimize
  1634. //  xlo = lower bound
  1635. //  xhi = upper bound
  1636. //  n = number of iterations (the bigger the more accurate)
  1637. fibmin = function(f, xlo, xhi, n) {
  1638.     local(a, b, x, y, ex, ey, k, lo, hi);
  1639.  
  1640.     lo = xlo;
  1641.     hi = xhi;
  1642.     k = n;
  1643.     for (k in k:2:-1)
  1644.     {
  1645.         a = fib(k - 2) / fib(k);
  1646.                 b = fib(k - 1) / fib(k);
  1647.                 x = lo + (hi - lo) * a;
  1648.                 y = lo + (hi - lo) * b;
  1649.                 ex = f(x);
  1650.                 ey = f(y);
  1651.                 if (ex >= ey)
  1652.                 {
  1653.                         lo = x;
  1654.                 else
  1655.                         hi = y;
  1656.                 }
  1657. //  printf("%d: (%g %g) %g %g %g %g\n",  k, a, b, lo, hi, ex, ey);
  1658.         }
  1659.     return (lo + hi) / 2;
  1660. };
  1661.  
  1662. //
  1663. // Simple example using above function to mimize f065. Answer is 0.65
  1664. //
  1665. x = fibmin(f065, 0, 1, 30); // printf("f(%g)=%g\n", x, f065(x));
  1666. if (abs(x - 0.65) > 1e-6)
  1667. {
  1668.     printf("x = %f\n", x);
  1669.     error("failed fibmin test");
  1670. }
  1671.  
  1672. printf("\tpassed fibmin test...\n");
  1673.  
  1674. //
  1675. //--------------------- Nasty Function Test ------------------------
  1676. //
  1677. printf("\tStarting Nasty Function Test...");
  1678. printf("\tthis will take awhile\n");
  1679. check = function( a, b, c, d, e, f, g, h ) {
  1680.   if ( a+b+c+d == e+f+g+h && ...
  1681.        a^2+b^2+c^2+d^2 == e^2+f^2+g^2+h^2 && ...
  1682.        a^3+b^3+c^3+d^3 == e^3+f^3+g^3+h^3 ) {
  1683.     return 1;
  1684.   else
  1685.     return 0;
  1686.   }
  1687. };
  1688.  
  1689. for(a in 8:10) {
  1690.   for(b in 7:(a-1)) {
  1691.     for(c in 6:(b-1)) {
  1692.       for(d in 5:(c-1)) {
  1693.         for(e in 4:(d-1)) {
  1694.           for(f in 3:(e-1)) {
  1695.             for(g in 2:(f-1)) {
  1696.               for(h in 1:(g-1)) {                  
  1697.               if(check( a, b, c, d,  e, f, g, h ) || ...
  1698.                      check( a, e, c, d,  b, f, g, h ) || ...
  1699.                      check( a, f, c, d,  e, b, g, h ) || ...
  1700.                      check( a, g, c, d,  e, f, b, h ) || ...
  1701.                      check( a, h, c, d,  e, f, g, b ) || ...
  1702.                      check( a, b, e, d,  c, f, g, h ) || ...
  1703.                      check( a, b, f, d,  e, c, g, h ) || ...
  1704.                      check( a, b, g, d,  e, f, c, h ) || ...
  1705.                      check( a, b, h, d,  e, f, g, c ) || ...
  1706.                      check( a, b, c, e,  d, f, g, h ) || ...
  1707.                      check( a, b, c, f,  e, d, g, h ) || ...
  1708.                      check( a, b, c, g,  e, f, d, h ) || ...
  1709.                      check( a, b, c, h,  e, f, g, d ) || ...
  1710.                      check( a, e, f, d,  b, c, g, h ) || ...
  1711.                      check( a, e, g, d,  b, f, c, h ) || ...
  1712.                      check( a, e, h, d,  b, f, g, c ) || ...
  1713.                      check( a, f, g, d,  e, b, c, h ) || ...
  1714.                      check( a, f, h, d,  e, b, g, c ) || ...
  1715.                      check( a, g, h, d,  e, f, b, c ) || ...
  1716.                      check( a, b, e, f,  c, d, g, h ) || ...
  1717.                      check( a, b, e, g,  c, f, d, h ) || ...
  1718.                      check( a, b, e, h,  c, f, g, d ) || ...
  1719.                      check( a, b, f, g,  e, c, d, h ) || ...
  1720.                      check( a, b, f, h,  e, c, g, d ) || ...
  1721.                      check( a, b, g, h,  e, f, c, d ) || ...
  1722.                      check( a, e, f, g,  e, f, g, h ) || ...
  1723.                      check( a, e, f, h,  e, f, g, h ) || ...
  1724.                      check( a, e, g, h,  e, f, g, h ) || ...
  1725.                      check( a, f, g, h,  e, f, g, h ) ) { cnt++; }
  1726.               }
  1727.             }
  1728.           }
  1729.         }
  1730.       }
  1731.     }
  1732.   }
  1733. }
  1734.  
  1735. if(1) {  // figure out the value of cnt, and check!
  1736.   printf("\tpassed nasty function test...\n");
  1737. else
  1738.   error();
  1739. }
  1740.  
  1741. //
  1742. //------------------ Test More Advanced Functions --------------------
  1743. //
  1744.  
  1745. printf( "\tStarting the lqr/ode test..." );
  1746. printf( "\tthis will take awhile\n" );
  1747.  
  1748. lqr = function( a, b, q, r ) {
  1749.     local( k, s,...
  1750.            m, n, mb, nb, mq, nq,...
  1751.            e, v, d );
  1752.  
  1753.     m = size(a)[1]; n = size(a)[2];
  1754.     mb = size(b)[1]; nb = size(b)[2];
  1755.     mq = size(q)[1]; nq = size(q)[2];
  1756.     
  1757.     if ( m != mq || n != nq ) {
  1758.     fprintf( "stderr", "A and Q must be the same size.\n" );
  1759.     quit
  1760.     }
  1761.     
  1762.     mr = size(r)[1]; nr = size(r)[2];
  1763.     if ( mr != nr || nb != mr ) {
  1764.     fprintf( "stderr", "B and R must be consistent.\n" );
  1765.     quit
  1766.     }
  1767.     
  1768.     nn = zeros( m, nb );
  1769.     
  1770.     // Start eigenvector decomposition by finding eigenvectors of Hamiltonian:
  1771.     
  1772.     e = eig( [ a, solve(r',b')'*b'; q, -a' ] );
  1773.     v = e.vec; d = e.val;
  1774.     
  1775.     index = sort( real( d ) ).ind;
  1776.     d = real( d[ index ] );
  1777.     
  1778.     if ( !( d[n] < 0 && d[n+1] > 0 ) ) {
  1779.     fprintf( "stderr", "Can't order eigenvalues.\n" );
  1780.     quit
  1781.     }
  1782.     
  1783.     chi = v[ 1:n; index[1:n] ];
  1784.     lambda = v[ (n+1):(2*n); index[1:n] ];
  1785.     s = -real(solve(chi',lambda')');
  1786.     k = solve( r, nn'+b'*s );
  1787.     
  1788.     return << k=k; s=s >>;
  1789.  
  1790. };
  1791.  
  1792. // Now run a little test problem.
  1793.  
  1794. k = 1; m = 1; c = .1;
  1795. a = [0     ,1    ,0    , 0;
  1796.     -k/m, -c/m,  k/m,  c/m;
  1797.      0,     0,    0,    1;
  1798.      k/m,  c/m, -k/m, -c/m ];
  1799. b = [ 0; 1/m; 0; 0 ];
  1800. qxx = diag( [0, 0, 100, 0] );
  1801. ruu = [1];
  1802. K = lqr( a, b, qxx, ruu ).k;
  1803.  
  1804. dot = function( t, x ) 
  1805. {
  1806.     return (a-b*K)*x + b*K*([1,0,1,0]');
  1807. };
  1808.  
  1809. x = ode ( dot, 0, 15, [0,0,0,0], .02, 1e-5, 1e-5 );
  1810.  
  1811. m = maxi( x[;2] );
  1812.  
  1813. if ( (abs( x[m;2] - 1.195 ) > 0.001)  || ...
  1814.      any (abs( x[x.nr;2,4] - 1 ) > 0.001) ) {
  1815.      printf( "\tfailed***\n" );
  1816. else
  1817.      printf( "\tpassed the lqr/ode test...\n" );
  1818. }
  1819.  
  1820. printf("Elapsed time = %i seconds\n", toc() );
  1821. "FINISHED TESTS"
  1822.